Quantum Monte-Carlo study of a two-species boson Hubbard model 
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We consider a two-species hard-core boson Hubbard model for a supersolid, where the two types 
of bosons represent vacancies and interstitials doped into a commensurate crystal. The on-site 
inter-species interaction may create bound states of vacancies and interstitials facihtating vacancy 
condensation at lower energies than in a single-species model, as suggested in an earlier mean field 
study. Here we carry out quantum Monte Carlo simulation to study possible supersolid phases of 
the model, corresponding to superfluid phases of the vacancies or interstitials. At low temperatures, 
we find three distinct superfiuid phases. The extent of the phases and the nature of the phase 
transitions are discussed in comparison to mean-field theory. 



I. INTRODUCTION 

A supersolid is a special type of solid with superfluid 
properties. It has a diagonal particle density long range 
order as in a usual crystal, and an off-diagonal long range 
order in particle density as in a superfluid. The simplest 
model for supersolid was proposed by Andreev and Lif- 
shitz in 1969.^ Their model was introduced to describe 
possible supersolid phase in Helium-4. In their model, 
vacancies or interstitials of solid Helium may exist in the 
ground state and condense due to the large quantum 
fluctuation of Helium atoms. The interaction between 
vacancy and interstital is neglected in their model. 

In this paper we study a two-species boson Hubbard 
model, which is an extension of the Andreev-Lifshitz 
model to include the interaction between vacancy and 
interstitial. This two-species model was recently intro- 
duced by Dai, Ma, and Zhang, ^ motivated by the obser- 
vation of non-classical rotational inertia moment in solid 
helium-4 reported by Kim and Chan.'^ They used a mean 
field theory to study the ground state of the model and 
the possibility of the supersolid phase. It was shown that 
the interaction of vacancies and interstitials may facili- 
tate a supersolid phase. In this paper we use quantum 
Monte Carlo (QMC) simulations to study the possible 
supersolid and the finite temperature phase transition 
in the two-species boson model. The simulations sup- 
port the qualitative conclusion obtained in the mean field 
theory that the vacancy-interstitial interaction may facil- 
itate supersolidity. Using QMC, we calculate the phase 
diagram, the superfluid densities of bosons, and the spe- 
cific heat of the system. The two-species boson model 
and our calculations may be useful to understand other 
boson problems such as bosons in optical lattices.'* 

Before we present the model and our results, we briefly 
summarize the current situation in study of supersolid 
Helium-4. Because of its light mass and its bosonic na- 
ture, solid helium-4 has been a natural candidate for pos- 
sible supersolid at low temperatures and high pressures. 
Theoretically, such a possibility was proposed by An- 
dreev and Lifshitz^ and by Chester Leggett further pre- 



dicted the non-classical rotational inertia moment of such 
a supersolid in a rotating experimenti^ii The interest of 
supersolid has been recently revived due to the observa- 
tion of non-classical rotational inertia in solid helium-4i^ 
By now, the non-classical inertia moment in solid helium 
has been confirmed by other groupsi^i^iiiS However, it re- 
mains controversial if the phenomenon is related to the 
supersolidity and if the supersolid phase is a bulk equilib- 
rium phenomenon.— On the theoretical side QMC 
simulations did not find a supersolid phase in Helium- 
^^ 14^15,16 Furthermore, the vacancies or interstitials in 
helium are shown to attract to each other and to tend 
to have phase-separation^ indicating that the Andreev- 
Lifshitz model may not describe solid helium. 



II. MODEL AND METHOD 

We consider a two-species boson Hubbard model in a 
cubic lattice with z = 6 nearest neighbors: 

3 

'^{taajaj + tbblbj + h.c.) 

where aj is an annihilation operator of boson a at lat- 
tice site j, representing a vacancy, and bj an annihila- 
tion operator of boson representing an interstitial, in a 
vacuum representing a defect-free insulating crystal of 
bosonic atoms, nj^a — and n^ f, = b'^bj are the num- 
ber operators for a- and 6-bosons, and Sa and £(, are site 
boson energies, respectively. We consider the interest- 
ing case > 0, and £(, > 0. We assume both vacancy 
and interstital are hard-core bosons, so that the allowed 
values for rij^a and n^.f, are either or 1. An exciton is 
described by the state with both a vacancy and an in- 
terstitial at the same lattice site rij^a — nj^h — 1. The 
couplings ta and ti, are the hopping integrals for boson a 
and b, respectively, and we assume ta and ti, to be posi- 
tive without loss of generality. U is the on-site attractive 
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interaction between a vacancy and an interstitial. Note 
that the attractive interaction between a vacancy and 
an interstitial reflects the strong short range repulsion 
between two nearby atoms when an interstitial atom is 
added into the lattice. 

Without interaction, at [7 = 0, the two-species model 
decouples into independent vacancy and interstitial mod- 
els. The ground state of the a-boson (vacancy) model 
is superfluid (a vacancy supersolid) if zta > Ca and an 
empty vacuum state (insulating solid) otherwise. Simi- 
larly, the ground state of the 6-boson (interstitial) model 
is superfluid (an interstitial supersolid) if zth > £{, and 
the empty vacuum state (an insulating solid) otherwise. 

The attractive inter-species boson interaction U cou- 
ples the two types of boson, and the problem can- 
not be solved analytically without approximation. This 
model was studied by using a mean field theory at zero 
temperature,^ and a special limiting case with Cb —>■ oo 
but a finite £b — U was investigated by a modified spin 
wave theory^ The main effect of the attractive term is 
to facilitate vacancy or/and interstital condensation due 
to excitons or bound states. 

In this paper we will use QMC methods to study the 
phase transition and superfluid properties of the model. 
We use a slightly extended version of the directed loop 
algorithmsiii^ of the ALPS projectfi^ In the stochastic 
series expansion (SSE) representation used by this algo- 
rithm the superfulid density can be measured through the 
fluctuations of the spatial winding numbers.-?. In three 
dimensions for a simple cubic lattice the relationship is: 

Ps,a = ^{WD, (2) 

where a — a,b refers to the type of boson, Wa is the spa- 
tial winding number of the bosons in one direction, L is 
the linear size of the cubic lattice and T the temperature. 
In addition we consider the correlated winding numbers 

W± = {{Wa ± Wb)^) = (W^) + (W^) + 2{WaWb) (3) 

and define 

P± = ^W±. (4) 

We performed simulations on lattices with up to 10'^ 
lattice sites. Larger lattices did not equilibrate using the 
directed loop algorithm due to the formation of bound 
states between a and b bosons. A two-worm algorithm 
such as developed for the one-dimensional case in Ref. 
[21I would be required to go to larger lattices, however we 
found that the sizes used here were sufficient to determine 
the nature of the phases. 

III. SYMMETRIC CASE 

We first consider the symmetric case with ea = £& = 1. 
The energy cost of an exciton, (both an a- and a 6-boson 
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FIG. 1: Finite temperature T phase diagram of model (1) in 
the symmetric case ta = tt = t, obtained by quantum Monte 
Carlo (squares with error bars), ea = tb = 1, U = 1.8. The 
dashed line is a linear fit to the data, separating supersolid 
phase (SS) from normal solid (NS). The mean field transition 
point at T = is indicated by a vertical arrow." Lines A - 
D indicate the parameters of cuts through the phase diagram 
that we will investigate in more detail in Figs. [2]and|4] 



on the same lattice site), is A = ea + eb ^ U . The sys- 
tem becomes a trivial exciton lattice if A < 0, which 
we will not discuss. For the symmetric case, we choose 
U — 1.8 in our simulations to examine correlation effects. 
In Fig. [Tl we summarize our result by showing the phase 
diagram in the parameter space of temperature T and 
boson hopping zt = zta = ztb- The simulations are car- 
ried out at temperatures ranging from 0.1 to 0.2, which 
allow us to estimate a zero temperature phase bound- 
ary at zt « 0.65, smaller than the mean field value of 
zt « 0.80, and smaller than that of the non-interacting 
case oi U = Q &i zt = 1. Hence, quantum fluctuations 
which are neglected in the mean field theory further favor 
the superfluid phase. 

We examine the temperature dependences in more de- 
tail along line A in the phase diagram of Fig. [T] Figure 
[2] shows the superfluid densities p± and the specific heat 
cy as functions of T along line A in Fig. [1] {zt — 0.88) 
at the temperature region from T — 0.03 to T = 0.4. 
As T decreases, p± rise abruptly below T — 0.13 and 
saturate to p± = 0.1, cy develops a clear peak around 
T = 0.13, and the peak becomes sharper as the size in- 
creases. Note that p-^ = p_ within our error bars, indi- 
cating that there are no correlations and the two types of 
bosons condense independently with the same superfluid 
density pa = Ph « P+/2. 

In Fig. [3l we show the the superfluid density for system 
sizes L = 4, 6, 8. Finite size scaling for a second order 
phase transition in the U{1) universality class implies 
that psL is a constant at the transition temperature Tc- 
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FIG. 2: Superfluid densities p± and specific heat cv as func- 
tions of tfie temperature T, at zta = ztb = 0.88 along line A 
in Fig. [U 




FIG. 3: Superfluid density psL as a function of T for the 
symmetry model at zt = 0.88, along line A in Fig[l] for L = 4, 
6 and 8. The three curves cross at one point, ihom which we 
estimate Tc = 0.129 ± 0.002. 



FIG. 4: Superfluid densities p± and specific heat as functions 
of zt in the symmetric model along line B in Fig. [1] 
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FIG. 5: PsL as functions of in the symmetric model at 
T — 0.1, 0.15, 0.20 corresponding to lines B, C and D shown 
in Fig. [1] For each temperature the curves for different system 
sizes cross at one point, consistent with a second order phase 
transition. 



The three curves in Fig. [3] indeed cross at a single point, 
from which we can estimate Tc — 0.129 ± 0.002. 

FinaUy we investigate the dependence on the hopping 
amplitude zt at various temperatures. In Fig. [4] we plot 
the superfluid densities and specific heat at T = 0.15 
(along line C). Superfluidity develops at around zt = 
0.92 as we can see from both p+ and cv- The superfluid 
density as functions of zt are plotted in Fig. [5]for different 
system sizes along lines B, C, and D. Each set of curves 
cross at one point, consistent with the expected scaling 
at a second order phase transition. 



IV. NON-SYMMETRIC CASE 

We now discuss the non-symmetric case, which is more 
interesting and possibly more relevant to physical sys- 
tems since there is a lack of vacancy-interstitial symme- 
try. In all the simulations reported for the non-symmetric 
model, we consider ea = l,eb = 4, and U — i. We use 
smaller values of U than in the mean-field work of Ref. 
12 since larger values of U cost too much CPU time. 

Our phase diagram in the parameter space of zta and 
zth at T = 0.15, is summarized in Fig. [Sfa). This tem- 
perature is low enough to observe the expected supersolid 
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FIG. 6: (a). Phase diagram obtained in QMC at T = 0.15 
for the non-symmetric case of model (1) with parameters U = 
4.0, Ca = 1 and eb = 4. The phase transition points marked 
by red circles are obtained from system sizes L = 3 and 4. 
Blue squares are obtained from systems with L = 4, 5 and 
6. Points A and B and the Lines C, D, E, F indicate cuts 
through the phase diagram that we study in more detail in 
Figs. [Tl llOl (b). Mean field ground state phase diagram for 
the same parameters, obtained using the mean field theory of 
Dai et al^ Dashed green lines are the transition lines of the 
noninteracting model with U = 0. 



phases. In addition to the normal solid, there are three 
supersolid phases: 

1. a vacancy superfluid-A phase [V-SF(A)] in which 
the a-bosons (vacancies) condense pa ^ and no 
6-bosons (interstitials) are present 

2. a vacancy superfluid-B phase [V-SF(B)] in which 
the 6-bosons condense pb ^ and ni^a = 1 (for 
T = 0). This is a vacancy superfluid above a back- 
ground of excitons. Vacancies move in an other- 
wise excitonic lattice, so it may be called vacancy 
superfluidi^ 

3. a vacancy and interstitial superfluid [VI-SF] phase 
in which both a- and 5-bosons condense: pa =/= 
and ph 0. 

The phase boundaries labeled by red circles in Fig.lHI^a) 
are obtained from simulations on systems with up to 
L = 4. Calculations on larger size systems in these pa- 
rameter region require much more computational effort, 
and are only carried out for four selected points, labeled 
by blue squares on the boundaries in the figure, represent- 
ing typical interesting cases of the three most interesting 
different phase-transitions in the parameter space. 

For comparison, we show in Fig. EUb) the result of 
mean field calculations for the same parameters consid- 
ered. Note that the QMC predicts a larger parameter 
space for the supersolid phases than the mean field the- 
ory, indicating again that the quantum fluctuation ne- 
glected in the mean field theory but included in the QMC 
is in favor of the supersolid phase. 
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FIG. 7: Superfluid density pb of fe-boson and specific heat as 
functions of T for non-symmetric case at point A in Fig. 6(a). 
Here the low-temperature phase is V-SF(B). 



In the remaining part of this section, we discuss the 
phase transition as a function of temperature and as a 
function of boson hopping integrals. 

To study the temperature dependence, we choose two 
typical points A [zta — 0.3 and ztf, — 3.75) and B 
{zta = 0.8 and ztf, = 3.5) in the parameter space as 
indicated in Fig. [6ja). In Fig. [71 we show the superfluid 
density and specific heat as functions of the temperature 
for the system at point A. As the temperature decreases, 
Pb starts to increase sharply at around T « 0.92, while 
Pa remains zero. This indicates that only 6— bosons con- 
dense. A scaling analysis of pbL gives Tc = 0.924 ±0.002. 

In Fig. [5] we show pa^b and cy for the system at the 
point B where there are two transitions. As the tempera- 
ture is lowered, the system first undergoes a transition at 
T = 0.814 ±0.002 from a normal-solid into the V-SF(B)- 
phase, in which the 6-bosons condense. As the temper- 
ature is further lowered, the system undergoes a second 
phase transition at T = 0.23 ± 0.01 within the super- 
solid state from the V-SF(B) phase to the VI-SF-phase 
where both the a- and 6— bosons condense. The criti- 
cal temperatures have again been estimated by a scaling 
analysis similar to the symmetric case. Note that be- 
low the lower transition point, pb further increases, due 
to the attractive interaction with the a-bosons which ef- 
fectively increase the chemical potential for the 6-bosons 
and hence their number. We have calculated p± and have 
found that p+ and p_ are almost the same, so that the 
correlations are very small. 

We now discuss the phase transitions along the lines C- 
F of Fig. 6(a) in more detail, at a temperature T — 0.15. 
There are three different phase-transitions: 

1. a transition between the insulating state and the 
V-SF(B) phase and along the line C in Fig. dja). 
For fixed zta = 0.2, we pass through a critical value 
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FIG. 8: Superfluid densities Pa and pi, and specific lieat as 
functions of T for non- symmetric case at point B indicated 
in Fig. [HI a). Here the low temperature phase is the VI-SF 
phase where both a-boson and &-boson condense. 
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FIG. 9: Superfluid densities pa and pb as functions of ztt for 
zta = 0.2 and zta = 0.7 along the lines C and D of Fig. |6]at 
T = 0.15. Note that pa = for zta = 0.2. 

where pb becomes non-zero while pa remains zero 
(see Fig. We estimate the critical value ztb = 
3.44 ± 0.02 using the scaling analysis. 

2. a transition between the insulating state and the 
VI-SF phase appears along the lines D and E in 
Fig. IDJa). For fixed zta = 0.7, the superfluid densi- 
ties for both a— and 6-bosons are zero at small val- 
ues of ztb, and become finite above a critical value, 
which is the same for the two types of bosons, as we 
can see from Fig. [HI The critical value of ztf, can 
be estimated using a scaling analysis for different 
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FIG. 10: Superfluid densities as functions of zta for non- 
symmetric case along line A in Fig. EJa) at T = 0.15. At 
this temperature, pb > 0, and the transition is between the 
two supersolid phases V-SF(B) and VI-SF. 



system sizes up to L = 6 which gives the critical 
values of ztb — 3.16 ± 0.02 from for the line D and 
ztb = 2.87 ± 0.02 for line E. 

3. a transition is between two supersolid phases along 
line F in Fig. [Ifa). Along this line ztb = 3.6 and 
Pb is always finite. As zta increases pa is zero up to 
a critical value zta — 0.55 ± 0.01. (see Fig. [TU|) . 



V. CONCLUSIONS 

Our quantum Monte Carlo simulations of a two-species 
bosonic Hubbard model of a supersolid show a phase di- 
agram qualitatively consistent with previous mean field 
results.^ The attractive interaction between a vacancy 
and interstitial may facilitate the superfluidity in a 
bosonic solid, even when single vacancies or interstitials 
are gapped. Quantum fluctuations which are ignored in 
the mean-field calculations stabilize the superfluid phase 
over a larger parameter regime. Unlike the modified spin 
wave calculations* which finds first order phase transi- 
tions at finite temperatures, the quantum Monte Carlo 
calculations show consistency with the expected scaling 
behavior at second order phase transitions in the U{1) 
universality class, for temperatures above T = 0.10. 
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